Magnetostatics solution (linear case)

%%capture
%run current_density.ipynb

##############################################################################
# Tree/Cotree gauging
##############################################################################
MESH = pde.mesh3.netgen(geoOCCmesh)
R = pde.tools.tree_cotree_gauge(MESH)

R

<27728x23597 sparse matrix of type '<class 'numpy.float64'>'
        with 23597 stored elements in Compressed Sparse Column format>

order = 1

D = pde.int.assemble3(MESH, order = order)

phix_Hcurl, phiy_Hcurl, phiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = order)
curlphix_Hcurl, curlphiy_Hcurl, curlphiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = order)

K_Hcurl = curlphix_Hcurl @ D @ curlphix_Hcurl.T +\
          curlphiy_Hcurl @ D @ curlphiy_Hcurl.T +\
          curlphiz_Hcurl @ D @ curlphiz_Hcurl.T

KR = R.T@K_Hcurl@R

r = jx_L2 @ D @ phix_Hcurl.T +\
    jy_L2 @ D @ phiy_Hcurl.T +\
    jz_L2 @ D @ phiz_Hcurl.T


cholKR = chol(KR.tocsc())
A = R @ cholKR.solve_A(R.T@r)

curlphix_Hcurl_P0, curlphiy_Hcurl_P0, curlphiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = 0)
Bx = curlphix_Hcurl_P0.T @ A
By = curlphiy_Hcurl_P0.T @ A
Bz = curlphiz_Hcurl_P0.T @ A

##############################################################################
# Storing to vtk
##############################################################################

grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_L2_Vector(grid,Bx,By,Bz,'grad_x')
pde.tools.vtklib.writeVTK(grid, 'magnetostatics_solution.vtu')

import pyvista as pv
mesh = pv.read('magnetostatics_solution.vtu')
mesh2 = pv.read('current_density.vtu')

mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,1])

p = pv.Plotter()
p.add_mesh(threshed, style='surface', color = "w", opacity=0.2, label=None)

threshed.set_active_vectors("grad_x")
arrows = mesh.glyph(scale="grad_x", orient=True, tolerance=0.03, factor=1000)
p.add_mesh(arrows, color="orange")


mesh2.set_active_vectors("J_L2")
arrows2 = mesh2.glyph(scale="J_L2", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows2, color="black")

p.camera_position = [(0, -500, 200),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend="html")